The nature of high-energy radiation damage in iron: Modeling results 
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Understanding and predicting a material's performance in response to high-energy radiation dam- 
age, as well as designing future materials to be used in intense radiation environments, requires the 
knowledge of the structure, morphology and amount of radiation-induced structural changes 1-5 . We 
report the results of molecular dynamics simulations of high-energy radiation damage in iron in the 
range 0.2-0.5 MeV. We analyze and quantify the nature of collision cascades both at the global and 
local scale. We find that the structure of high-energy collision cascades becomes increasingly con- 
tinuous as opposed to showing sub-cascade branching reported previously. At the local length scale, 
we find large defect clusters and novel small vacancy and interstitial clusters. These features form 
the basis for physical models aimed at understanding the effects of high energy radiation damage 
in structural materials. 



PACS numbers: 

I. INTRODUCTION 

Radiation effects are common in nature. Their sources 
vary from cosmic radiation to decay of isotopes in ter- 
restrial rocks. A large variety of radiation sources are 
also created and used in science and technology. This 
includes energy generation in existing nuclear power sta- 
tions, where kinetic energy of fission products is con- 
verted into heat and electricity. When feasible, future fu- 
sion reactors will harvest the energy from thermonuclear 
reactions. In these applications, the energy of emitted 
particles has a two- fold effect: on one hand, this energy 
is converted into useful energy, by heating the material; 
on the other hand, this energy damages the material and 
degrades the properties important for the operation, in- 
cluding mechanical, thermal, transport and other proper- 
ties. This is currently a central issue for fusion reactors, 
where the ability of metal structural components to with- 
stand very high neutron fluxes is intensely discussed 1 3 . 
Another example is the damage to nuclear reactor ma- 
terials coming from fission products. In addition, the 
nuclear industry faces yet another problem, that of ra- 
diation damage to materials to be used to encapsulate 
long-lived radioactive waste 4,6,7 . 

A heavy energetic particle displaces atoms on its path 
which, in turn, displace other atoms in the system. A 
collection of these atoms is often referred to as a "colli- 
sion cascade" 8 14 . A typical collision cascade created by 
a heavy 100 keV particle propagates and relaxes on the 
order of picoseconds and spans nanometers. The result- 



ing structural damage, in the form of amorphous pock- 
ets or point defects and their clusters, ultimately defines 
to what extent materials' mechanical, thermal and other 
properties are altered. For example, radiation-induced 
defects can reduce materials' thermal conductivity and 
therefore result in inefficient energy transfer in both fu- 
sion and nuclear reactors, heat localization and other un- 
wanted effects. 

Understanding and predicting these and other effects, 
as well as designing future materials to be used in in- 
tense radiation environments, requires developing physi- 
cal models. These, in turn, are based on the knowledge of 
what high-energy radiation damage is in terms of struc- 
ture and morphology. 

Molecular dynamics (MD) simulations have been an 
important method for studying radiation damage in ma- 
terials because they give access to the small time and 
length scales of the collision cascades, and give a detailed 
picture of the damage at the atomistic scale. Previous 
MD simulations have provided important insights into 
the radiation damage process 8 15 . However, due to sys- 
tem size limitations in MD simulations, the reported re- 
sults were limited to energies of about 100 keV. Recently, 
the number of displaced atoms was reported for 200 keV 
cascade 16 . 

On the other hand, knock-on energies are larger in sev- 
eral important applications. When impacted on 14 MeV 
neutrons, iron knock-on atoms in fusion reactors reach 
the energy of up to 1 MeV 2,9,16 with an average energy 
of about 0.5 MeV 15 . In fission nuclear reactions, the 
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fission product energies are on the order of 50 to 100 
MeV, transferring high energy to the surrounding mate- 
rial. The need to simulate realistic energy cascades has 
been particularly emphasized, with a view that extrapo- 
lation of low-energy results may not account for some im- 
portant features of higher-energy radiation process which 
can contain novel qualitative features. More generally, 
the need to simulate length and energy scales that are 
relevant and appropriate to a particular physical process 
has been recognized and reiterated 5 . 

In this paper, we study the radiation damage process 
due to high-energy Fe knock-on atoms of 0.2-0.5 MeV 
energy. We focus on high-energy radiation damage in a- 
iron, the main structural material in fusion and future 
fission reactors. We analyze and quantify the nature of 
collision cascades both at the global and local scale. We 
find that high-energy collision cascades may propagate 
and relax as increasingly continuous damage structures 
as opposed to showing sub-cascade branching as assumed 
previously. At the local length scale, we find large clus- 
ters and new defect structures. 



II. METHODS 

We have used the DL_POLY program, a general- 
purpose package designed for large-scale simulations 17 ' 18 . 
We have simulated systems with 100-500 million of 
atoms, and run MD simulations on 20000-60000 paral- 
lel processors of the HECToR National Supercomputing 
Service 19 . MD simulations were performed at a constant 
energy and volume ensemble (except for the boundary 
layer atoms, see below) with an initial temperature set 
to 300 K, which was preceded by equilibration runs in a 
constant pressure ensemble at 300 K. Periodic boundary 
conditions were imposed in all directions. 

We have implemented several features to handle radi- 
ation damage simulations. First, we used a variable time 
step to account for faster atomic motion at the begin- 
ning of the cascade development and its gradual slowing 
down at later stages. Second, the MD box boundary 
layer of thickness of about 10 A was connected to a con- 
stant temperature thermostat at 300 K to emulate the 
effect of energy dissipation into the sample. Third, we 
have accounted for the electronic energy losses (EEL), 
particularly important at high energies. EEL is a com- 
plicated process that involves a wide range of effects af- 
fecting damage production and annealing 2,20 22 . Tak- 
ing EEL into account in MD simulations involves mod- 
els that include slowing-down of atoms due to energy 
transfer to electron excitation processes as well as feed- 
ing this energy back to the system as the phonon energy. 
The implementation of this dual energy exchange mech- 
anism in MD simulations, based on the two-temperature 
approach 21,23 ' 42 , is in progress. In this work, we model 
electronic energy loss as a friction term added to the 
equations of motion. The characteristic energy loss re- 
laxation time (taken here as r es = 1.0 ps), is obtained by 



relating the stopping strength (A = 0.1093 eV 1 / 2 A~ 1 ) 26 
in the low-velocity regime via Lindhard's model to the 
rate of energy loss for a single atom 23,43 . Such electronic 
stopping would only be effective above a certain thresh- 
old, where the atoms would have sufficient energy to scat- 
ter inelastically. We use a cut-off kinetic energy value 
(8.6 eV) corresponding to twice of the cohesive energy 44 , 
however a number of other threshold values have been 
proposed 45 48 . 

For a-Fe, we have used a many-body embedded-atom 
potential 24 , optimized for better reproduction of several 
important properties of a-Fe including the energetics of 
point defects and their clusters ("M07" from Ref. 25 ). At 
distances shorter than 1 A, interatomic potentials were 
joined to short-range repulsive ZBL potentials 26 . 

To analyze the collision cascade, we have employed two 
methods. First, an atom is identified as "displaced" if it 
moves more than distance d from its initial position. The 
number of displaced atoms, iVdi sp , quantifies the overall 
amount of introduced damage. Some of this damage re- 
covers back to crystal. To account for this effect, we 
employed the second method in which an atom is identi- 
fied as a "defect" if it is further than distance d from any 
of the crystalline position of the original lattice. Then, 
the number of defects TV^ef quantifies both damage pro- 
duction and its recovery. We note that A^isp and TVdef 
depend on d (we use d = 0.75 A as a convenient measure). 
However, the trends discussed in Sec. Ill, including the 
two regimes of cascade relaxation as well as dynamics of 
defects recovery are not sensitive to d provided it is in the 
sensible range of distances (e.g., too small d < 0.1 — 0.2 
A will be affected by usual thermal fluctuations whereas 
d> 1 - 1.5 A may not identify defect atoms). 

Vacancies or self-interstitial atoms (SI A) are defined 
to belong to the same defect group (cluster) if within 
2 nearest-neighbour distance (plus a 20% perturbation). 
Second nearest-neighbour (nn) is a common clustering 
criterion for SIAs 16 ' 36 , however the criteria for vacancy 
clusters vary significantly (from 1 to 4nn) across the 
literature 38 . When identifying cluster size in the Sec. Ill 
the net defect number (the difference between the num- 
ber of SIAs and the number of vacancies) is reported. 



III. RESULTS AND DISCUSSION 

We discuss the main features of high-energy collision 
cascades. To account for potentially different collision 
cascades due to different knock-on directions, we have 
simulated 4 different directions for each energy, avoiding 
low-density directions and associated channeling. 

TVdisp and TVdef are shown in Fig. 1 for 0.2 and 0.5 MeV 
cascades simulated in different knock-on directions. We 
observe two types of cascade relaxation in Fig. 1. The 
first type is related to the large peak of N^ sp « 10 6 at 
short times of about 1-2 ps. This peak relaxes during 
about 10 ps. This peak is often discussed as the "ther- 
mal spike" 28,29 related to melting inside the collision cas- 
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FIG. 1: iVdisp and iV def from 0.2 MeV (top) and 0.5 MeV 
(bottom) knock-on atoms for four cascades. 



cade and associated swelling that causes the temporary 
increase of A^isp- At the atomistic level, the increase of 
^Vdisp can be understood on the basis of anharmonicity 
of interatomic interactions: large-scale atomic motion in- 
side the cascade causes the increase of interatomic sepa- 
rations due to anharmonicity. This results in the outward 
pressure from the cascade on the surrounding lattice and 
lattice expansion. This elastic deformation lasts several 
ps, equal to several periods of atomic vibrations during 
which the energy is dissipated to the lattice, and gives 
rise to the peak of N^ sp . Notably, the elastic deforma- 
tion is reversible irrespective of whether it is followed by 
the recovery of the structural damage discussed below. 
In Fig. 1, A^isp, averaged over all knock-on directions at 
the end of simulation (corresponding to the flat lines in 
Fig. 1) is about 67,000 and 111,000 atoms for 0.2 MeV 
and 0.5 MeV cascades, respectively. 

Figures 2-3 show the snapshots of the collisions cas- 
cades at three different stages of cascade development 
that include the intermediate stage corresponding to the 
large peak of A^isp in Fig. 1. Consistent with Fig. 1, we 
observe a significantly larger number of atoms involved 
in large displacements at intermediate times as compared 
to the final relaxed state. 

The second type of cascade relaxation is related to the 
dynamics of A^ef- At short ps times, the large peak of 
TVdef is of the same origin as that seen for A^isp- However, 
dynamics of A^ef also reflects the recovery of structural 
damage. This recovery proceeds by the diffusion and re- 
combination processes during which atoms settle at the 
newly found crystalline positions. This process lasts up 
to 20 ps, significantly longer than relaxation time of the 
first elastic relaxation process (see Fig. 1). As a result 
of this relaxation, A^ef, averaged over all simulated di- 
rections at the end of simulation (corresponding to the 
flat lines in Fig. 1), is about 1,800 and 2,800 atoms, re- 
spectively, corresponding to approximately 97% recovery 




FIG. 2: Displaced (top) and defect (bottom) atoms in a rep- 
resentative 0.2 MeV collision cascade. The knock-on atom 
moves from the top left to the bottom right corner. The three 
frames for each type of atoms correspond to 0.3 ps, 3 ps and 80 
ps, respectively. Cascade size (maximal separation between 
any two atoms in the cascade) is 560 A. Vacancies (inter- 
stitials) are represented in purple (green); we used Atomeye 
software 27 to visualize cascade evolution. 






FIG. 3: (a) and (b) show two representative 0.5 MeV cas- 
cades. The knock-on atom moves from the top left to the 
bottom right corner. In (a) we show both displaced (top) 
and defect (bottom) atoms at 0.2, 1.5 ps and 100 ps. In (b) 
we show the defect atoms only at 0.3 ps, 2 ps and 100 ps. 
Cascade size in (a) and (b) is 950 A and 1300 A, respectively. 
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rate as compared to A^isp- Such a high recovery rate is in 
interesting resemblance to some of the resistant ceramic 
materials, but in contrast to others 30,31 . 

Our simulations provide an important insight about 
the structure and morphology of high-energy cascades. 
The existing view of the high-energy cascade is that it 
branches out to smaller sub-cascades and "pockets" of 
damage that are well separated from each other. This 
takes places over a certain energy threshold, even though 
this threshold was not firmly established 16,32,33 . This 
picture originated as a result of using binary-collision 
simulations in combination with MD simulations of low- 
energy events. Although involving approximations in- 
herent in binary-collision simulations and extrapolations 
of low-energy MD simulations, this picture would offer a 
great degree of reduction and simplification: in analyz- 
ing the results and consequences of high-energy damage, 
only small sub-cascades need to be considered. 

In Figures 2-3, we observe a qualitatively different 
picture. Cascades branching is visibly reduced as com- 
pared to low-energy events, in that we do not find well- 
separated sub-cascades. Some cascade branching is seen 
in the first 0.5 MeV cascade shown in Fig. 3a only and, 
importantly, during an intermediate stage of cascade de- 
velopment only. On the other hand, the final cascade 
morphology is described by a rather continuous distri- 
bution of the damage across about 1000-1300 A where 
no distinctly separated sub-cascades can be identified. 
Common to all collision cascades we have simulated, this 
picture is particularly visible for defect atoms in the final 
state of the cascades shown in Figs. 2-3. 

Qualitatively, reduced cascade branching and the 
emergence of a more continuous damage distribution can 
be understood as follows. For a scattered atom to move 
far enough from its initial position and form a spatially 
separated sub-cascade (i.e., branch out) requires a large 
value of energy transferred to it by the incident atom. In 
the absence of inelastic losses, the transferred energy, T, 
is T = |T m (l — cos (0)), where T m is the maximal trans- 
ferred energy and <fi is the scattering angle 34 . For large 
energy of the incident atom, E, (j) decreases as <j) oc -^ 34 . 
We therefore find that for large E and small </>, T de- 
creases as T oc (j) 2 ex -^2. Large E and, consequently, 
small T, results in scattered atoms forming the damaged 
region close to the trajectory of the initial knock-on atom 
and, therefore, promote a continuous structure of the re- 
sulting damage. This is consistent with our current find- 
ings. 

We note that as the incident atom slows down, T in- 
creases, leading to sub-cascade branching at the end of 
the atom trajectory. However, larger E results in the in- 
crease of the relative fraction of continuous damage over 
the fraction of branched cascades. 

Our finding is important in the context of the long- 
term evolution of radiation damage. Indeed, a recent 
kinetic Monte Carlo study 36 has shown that very large 
defect clusters can have a major effect on the long-term 
damage development. 



The discussion has so far concentrated on the large- 
scale cascade structure and morphology. We now briefly 
discuss defect structures at the local level. The size and 
structure of the defect clusters created by the cascades 
were analyzed and the results are summarized in Table I. 
The simulations confirm that the normalized fraction of 
Frenkel pairs (FPs) of 0.3-0.4 is roughly constant for cas- 
cades over 0.1 MeV 16,38 . The number of surviving FPs 
follows the trend from lower cascade energies (Fig. 4). 
The fraction of surviving interstitials grouped into clus- 
ters was found to be 0.58(3) and 0.52(3) for cascade en- 
ergies of 0.2 and 0.5 MeV, respectively. This is consistent 
with the results for 50-100 keV cascades 16,36 hinting at a 
possibility that the clustering fraction may reach a max- 
imum at ~ 100 keV. A similar trend was observed for 
vacancy clustering fractions of 0.33(3) and 0.35(1), re- 
spectively. 
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FIG. 4: Dependence of the number of surviving FPs on 
the cascade energy. A comparison between Finnis- Sinclair 
potential runs started at 100 K (red line 16 ) with the current 
work (blue). The fits correspond to a single-cascade and sub- 
cascades (branching) regimes (details in 16 ). 

The majority of interstitial clusters were found to be 
glissile <111> crowdion clusters. The largest defect 
structure was a composite 89-interstitial cluster, formed 
from a merger of a set of <111> and <100> crowdions 
(Fig. 5). Owing to its complex morphology, it will be 
immobile. This interstitial cluster is quite large, yet con- 
sistent with the data reported for lower energies 35,36 . 

Large vacancy clusters were also observed with the 
54-vacancy cluster being the largest one. Several large 
vacancy clusters formed <100> and <111> dislocation 
loop-like configurations [Fig. 6(a)]. A cross section of an 
exemplar vacancy cluster is shown [Fig. 6(b)], to empha- 
size its dislocation-like nature. The smaller vacancy clus- 
ters revealed a rich variety of structures, such as hexag- 
onal vacancy clusters with interstitial rings surrounding 
a central vacancy [Fig. 7(a)]. 

We also observed a wide range of sessile interstitial 
clusters. Some of these could be clearly identified as be- 
ing related to the C15 meta-stable phase discussed in 41 , 
and many were joined to crowdions or crowdion clusters 
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1450 (220) 


0.29 (0.04) 


150 (14) 
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36 (7) 
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47 
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TABLE I: The number of Frenkel pairs (FP), (N-p), and defect statistics for 0.2 MeV and 0.5 MeV cascade simulations in 
a-iron. The value in brackets shows standard error in the mean over 4 constituent runs for each simulation. Largest clusters 
are determined by net defect count. NRT fraction is the normalized number of FP 37 . 




FIG. 5: The largest cluster consisting of 89 intersitials. It is 
mainly composed of a set of <100> crowdions (selected region 
in (a)) and a fraction of normally glissile <111> crowdions 
(region highlighted in (b)). Such cluster morphology blocks 
the motion of crowdions and results in an overall sessile clus- 
ter; similar effect of immobilization of a cluster by another 
defect was observed in 39 . Interstitials (vacancies) are shown 
in silver (blue). We used VMD package for visualization of 
local defect structures 40 . 




which we will need to be included in physical models 
aimed at understanding and predicting the effects of radi- 
ation damage on the mechanical, thermal and transport 
properties of structural materials. The reported dam- 
age structures such as the increased continuous morphol- 
ogy of high-energy collision cascades will form a starting 
point for long-timescale models in order to understand 
and predict the effects of radiation damage. Large defect 
structures reported here, including novel vacancy and in- 
terstitial clusters, will be important for understanding of 
the interaction of these clusters with transmutation gases 
and nucleation of helium bubbles. 




FIG. 7: (a) The (111) projection of a small 9-vacancy cluster; 

(b) a C15 phase tetra-interstitial with a <111> crowdion at- 
tached. The vacancies are omitted from the figure for clarity; 

(c) a hexagonal di-interstitial with a split interstitial attached. 



FIG. 6: (a) A (111) projection of a 39-vacancy cluster; (b) 
a (001) projection of this cluster. The large spheres show the 
central vacancy for selected constituent 'vacancy crowdions', 
thus emphasizing the dislocation nature of the cluster. 

[Fig. 7(b)]. Smaller ring-like structures were also ob- 
served [Fig. 7(c)] in which 6 atoms shared 4 neighbouring 
lattice sites. 

IV. CONCLUSION 

In summary, we have found novel structural features 
of radiation damage in iron on both large and local scale 
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